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The liquid regime equation of state of silicon dioxide Si02 is calculated via quantum molecular 
dynamics in the density range 5 to 15 g/cc and with temperatures from 0.5 to 100 eV, including the 
Q-quartz and stishovite phase Hugoniot curves. Below 8 eV calculations are based on Kohn-Sham 
density functional theory (DFT), above 8 eV a new orbital-free DFT formulation, presented here, 
based on matching Kohn-Sham DFT calculations is employed. Recent experimental shock data 
is found to be in very good agreement with the current results. Finally both experimental and 
simulation data are used in constructing a new liquid regime equation of state table for Si 02 . 


Access to the accurate equation of state (EOS) of ma¬ 
terials over wide ranges in density and temperature, and 
in particular for matter in extreme conditions of signifi¬ 
cantly elevated temperature and density with respect to 
ambient conditions, is important in various arenas. Some 
systems of interest include dense astrophysical plasmas 
as exist in the interiors of giant planets, as well as warm 
dense matter, which is increasingly studied in high en¬ 
ergy density laboratory experiments, and also in devel¬ 
opment of inertial confinement fusion [ij. A particular 
case, presented here is that of quartz (Si 02 ), a ubiqui¬ 
tous mineral in Earth’s composition, it is also constituent 
in exoplanet modeling at higher temperatures and pres¬ 
sures than at inner Earth conditions ^. This material is 
also important as a window material in high compression 
shock experiments, typically in the quartz liquid regime 
above 100 GPa and 5000 K [^, |^. Hence, an accurate 
EOS for quartz is critical for determining properties of 
other materials through shock experiments. 

In general these extreme conditions represent a signifi¬ 
cant and current challenge of high energy density physics. 
Experimental results remain sparse, and from a theo¬ 
retical standpoint the ions exhibit moderate to strong 
coupling while the electrons require quantum treatment 
0. This then necessitates numeric simulations for which 
quantum molecular dynamics (QMD) based on Kohn- 
Sham density functional theory (DET) has emerged as 
the state of the art. Kohn-Sham DFT, however, becomes 
computationally prohibitive with increasing temperature 
as the number of required orbitals increases with tem¬ 
perature and in general the method scales as the cube 
of the number of orbitals. Though material and den¬ 
sity dependent, Kohn-Sham calculations are, thus, gen¬ 
erally viable below 10 eV or so. In this work we develop 
an orbital-free DFT formulation to combine with Kohn- 
Sham results and extend QMD simulations to very high 
temperatures. 

In QMD, the ions are treated classically and moved 
according to Newton’s equations, where the force on each 
ion is found from the Coulomb repulsion between all ions, 
and from the neutralizing electron charge density. The 
electron density, n, is found at each ionic configuration by 


DFT. This is done by minimizing the free energy, which 
is given by the density functional [1] 


F[n] = Fs[n] + [n] -|- F,cc[n] + Fei[n\ 


( 1 ) 


where F^ is the non-interacting free energy comprised of 
both kinetic and entropic parts, Fh is the Hartree energy 
or direct Coulomb interaction between the electrons, F^i 
is the electron-ion Coulomb interaction, and F^c is de¬ 
fined as the remainder of the total free energy, which in¬ 
cludes the quantum mechanical exchange and correlation 
as well as the excess kinetic and entropic terms. Of the 
contributions neither Fs nor F^c have exact formulations 
in terms of the density alone. Civen the same orbital- 
free Fxc approximation, the only difference in approach 
of orbital-free DFT from Kohn-Sham DFT is that the 
non-interacting free energy, F^, is found from an approxi¬ 
mate density functional instead of being exactly obtained 
through the calculation of single particle orbitals. 

In recent years the orbital-free approach at finite tem¬ 
perature has gained attention, with most results being 
for hot dense systems where the Thomas-Fermi approx¬ 
imation is employed for Fs Q- Various works have of¬ 
fered density gradient corrections to Thomas-Fermi that 
improves results moderately iSIl. None of these func¬ 
tionals, though, have reached the accuracy of Kohn-Sham 
across temperature regimes. A recent nonlocal functional 
has been shown to be hig hly accurate across temperature 
and density regimes [l2j |. but is requisite on Kohn-Sham 
derived pseudopotentials, which may not be as transfer¬ 
able. Subsequently an accurate and general orbital-free 
funcitonal remains elusive. 

Previous works [H, [l3| have attempted to connect 
high temperature Thomas-Fermi calculations with low 
temperature Kohn-Sham results, for low atomic num¬ 
ber systems. In this work we develop and implement 
a method to obtain a simple and accurate Fs beyond 
Thomas-Fermi, applicable to a wide range of materials 
and densities, for warm to hot systems, which smoothly 
extends Kohn-Sham results beyond the lower tempera¬ 
ture region that is currently computationally accessible. 
Then in an application to Si02 a new liquid phase EOS 
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is developed utilizing the current QMD calculations and 
compared with experimental results. 

We first introduce an orbital-free Fg of the following 
form 


= 

Ftf N + ^Fy\Y [n] 

(2) 

Ftf [n] = J 

^ F^^°{n{v),T)/Vdv , 

(3) 

Fvw N = J 

f|Vn(r)p^_. 

8n(r) 

(4) 


Here Ftf is the Thomas-Fermi approximation which 
takes for the noninteracting free energy that of the ho¬ 
mogeneous electron gas, with density equal to 

that of the local density and at the system temperature 
T ; Fyw is the von-Weizsacker gradient correction which 
has no explicit temperature dependence. The von Weiz- 
sacker coefficient A, is material and density dependent, 
and is to be determined through matching conditions 
between orbital-free and Kohn-Sham calculations at a 
specific density. In a more general way one could write 
A (no), with no being the average electron density, here 
though, A is treated as a constant and in that way Eq. 
0 represents a best approximation at a given density 
to the more general Fg, which is a single functional over 
varying densities. 

The procedure to determine A is straight forward. First 
we evaluate the pressure along an isochore using Kohn- 
Sham MD, until the temperature is high enough that 
the calculations becomes intractable (which of course de¬ 
pends on computational resources) but should be at least 
5 to 10 eV. Below which the orbital-free approach will be¬ 
come inaccurate due to issues such as molecular bonding. 
Then at a given match temperature a few orbital-free cal¬ 
culations are performed with initial guesses for A and the 
A is determined which reproduces the Kohn-Sham data 
in pressure, A is then fixed for that density. Figure [T] 
shows results of the matching method over a range of 
densities, a strong argument is made for the approach in 
that even though the match for pressure is made at the 
single point T = 6 eV, the change with temperature, or 
slope in Fig. [U is in near exact agreement between the 
two DFT methods. Also at very high temperature ^100 
eV, our results come into agreement with the Thomas- 
Fermi based MD, which is correct in the high-T limit. 

The details of the calculations are as follows. The 
Kohn-Sham Si02 calculations were performed in the liq¬ 
uid regime above 5.5 g/cc and 0.5 eV using the Quantum- 
Espresso program [l^. We included 72 total atoms in 
the calculations and all were performed with the F-point 
only. Some calculations were performed with a 2x2x2 k- 
point grid, without noticeable differences to the Gamma- 
point calculations. PAW pseudopotentials [1^ were em¬ 



FIG. 1: Pressure results of Kohn-Sham (KS) MD and orbital- 
free (OF) MD near the match temperature of 6 eV, for iso¬ 
chores of 11, 10, 9, 8, 6.9, 6.43, and 5.57 g/cc, A ranges from 
0.18-0.22. 


ployed, along with a temperature dependent implementa¬ 
tion of the AM05 exchange-correlation functional l3, 3 ; 
and a plane wave energy cutoff of 50 Ry was used. The 
molecular dynamics were performed at constant temper¬ 
ature using the Anderson thermostat. In order to calcu¬ 
late the Hugoniot, described below, the initial crystalline 
states of a-quartz and stishovite were calculated at the 
experimentally prescribed densities of 2.65 and 4.29 g/cc 
at room temperature, using the primitive cell of each 
phase with fully converged fc-point grids. 

In the orbital-free case, the same number of atoms were 
included as for the Kohn-Sham case, and the electron 
density was optimized on a regular 64^ numeric grid. Lo¬ 
cal pseudopotentials were generated for each density and 
temperature, according to the prescription given in Ref. 
0, with a cutoff radius of 0.6 time the Wigner-Seitz ra¬ 
dius. Since the orbital-free calculations were performed 
only above 5 eV, only the temperature dependent local 
density approximation (LDA) exchange-correlation [3 
was used for simplicity. The orbital-free molecular dy¬ 
namics were completed in the isokinetic ensemble (3 - 

The key application of the above QMD calculations 
is in construction and validation of a far reaching EOS 
valid at arbitrary densities and temperatures. While ex¬ 
perimental data serves as the traditional constraining in¬ 
put for EOS construction, it is particularly lacking in the 
warm dense matter regime. Further this regime falls in 
the region of interpolation between the low and high tem¬ 
perature models upon which the EOS is built, precisely 
where constraining input is most needed. The current 
QMD calculations then provide that input. 

The overall EOS model itself is a summation of three 
independent terms for the total Helmholtz free energy, 
given in terms of the material density p and the temper¬ 
ature T, 


Ftot — Fo{p) + Fi{p, T) + Fe{p, T) . (5) 
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The contributions here are the cold curve, Fq, the ion 
thermal contribution Fi, and the electron thermal con¬ 
tribution Fg. First of all the cold curve, representing 
zero temperature electrons and ions, is constructed from 
a modified Lennard-Jones model in the expanded region 
and Thomas-Fermi-Dirac theory at very high pressures, 
with a Birch-Murnaghan model in the interim. Next 
the ion thermal portion, due to the thermal motion of 
the ions, is built from a Debye model for temperatures 
less than the melt temperature, and interpolating to an 
ideal gas at high temperatures, leaving a poorly described 
liquid-like warm dense matter state. Finally the electron 
thermal portion, which comprises the thermal excitation 
of electrons above the ground state, is evaluated via the 
Thomas-Fermi average atom model. 

In order to effect an accurate EOS, the various param¬ 
eters of the constituent models need to be constrained by 
input data from experiment or further theoretical calcu¬ 
lations. 

In the case of Si02 there is shock data providing the 
pressure-density curve of the Hugoniot for a-quartz Q 
and more recently for stishovite as well as some recent 
shock release data for a-quartz which provides some 
off-Hugoniot data. This data alone is inadequate for a 
full ranging EOS, but when combined with wide ranging 
QMD data a more complete view is established, and an 
accurate EOS may be determined. 

The Hugoniot itself is calculated via the Rankine- 
Hugoniot jump conditions [2l|, which relates the equi¬ 
librium pre-shock and post-shock states by consideration 
of conservation laws. 


E-Eo = {P + Po){Vo-V)/2. (6) 

Here E, P and V are the internal energy, pressure and 
volume respectively, and the 0 subscript denotes the ini¬ 
tial state. For a-quartz and stishovite the initial condi¬ 
tions are densities of 2.65 g/cc and 4.29 g/cc respectively, 
and at ambient pressure and temperature. The initial 
internal energy is calculated via Kohn-Sham DFT and 
then using Kohn-Sham MD, an isochore or isotherm is 
followed until conditions are such that Eq. m is satis¬ 
fied. However when extending to orbital-free MD there is 
a shift in the internal energy due to change in the pseu¬ 
dopotential and the functional. As part of the matching 
scheme this shift is found at the match point temperature 
and density, so that the initial state energy calculated by 
Kohn-Sham DFT can be used. Additionally since there 
is a different pseudopotential used at each temperature 
and density with the orbital-free DFT, a second shift in 
energy must be accounted for, which is found by perform¬ 
ing two orbital-free calculations at the same temperature 
and density but with different pseudopotentials. 

The resulting Hugoniots for a-quartz and stishovite 
are shown in Fig. [2j The very good agreement between 
experiment and QMD, provides validity to the QMD re¬ 
sults for extension to regions where there is no experi- 
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FIG. 2: Comparison of experimental results, Refs. 13-11, with 
current QMD and EOS shock Hugoniot results. Upper panel 
shows Q-quartz, middle panel shows stishovite, and the lower 
panel shows a-quartz on the left curve and stishovite to the 
right. 


mental data, such as the high pressure Hugoniot exten¬ 
sions shown in the lower panel of Fig. [5] A benefit of 
the QMD calculation is that the temperature is also cal¬ 
culated along the Hugoniot. This is not always available 
from the shock data, as it is not here with the a-quartz 
data. Along the a-quartz Hugoniot we find the temper¬ 
ature increases from 0.86 eV at 5.57 g/cc to 6.73 eV at 
7.62 g/cc to 23.34 eV at 10.00 g/cc. 

Away from the Hugoniot the broad range of QMD re¬ 
sults can be compared directly with the EOS. A subset of 
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FIG. 3: Comparison of pressures between the QMD results 
and the current EOS for liquid Si02 isotherms. From top to 
bottom temperatures are 100, 50, 20, 10 ,5 eV in the upper 
panel and 6.89, 5.17, 3.45, 1.72, 0.86 eV in the lower panel. 
The a-quartz (dashed) and stishovite (dotted) Hugoniots cal¬ 
culated from the EOS are shown for reference. 

QMD calculations is plotted in Fig. |3] showing the range 
of temperatures for 0.86-100 eV and densities from 5-15 
g/cc. The higher temperature orbital-free results were 
instrumental in constraining the ion thermal portion of 
the EOS through the higher temperature liquid phase 
approaching the Thomas-Fermi-Dirac limit. 

Finally we show the agreement of the resulting EOS 
with the shock release data in Fig. ID Here the adiabats 
are calculated within the EOS which pass through the 
experimental release points along the Hugoniot. Good 
agreement is shown with experiment for the release in 
each of the three materials, which also with the isotherm 
data of Fig. |3] demonstrate high accuracy for the EOS 
away from the Hugoniot. 

In summary we have developed a robust and pre¬ 
scriptive method for accurately extending Kohn-Sham 
DFT based molecular dynamics simulations by orbital- 
free DFT simulations, for equation of state calculations 
including calculation of the Hugoniot to very high tem¬ 
peratures. This is important as the Kohn-Sham MD 
calculations have been shown to accurately characterize 
both solid and liquid systems and have become a gold 
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FIG. 4: Good agreement is seen between the experimental Z- 
machine shock release data. Ref. 0, and the EOS calculated 
adiabats. 


standard. Yet due to temperature scaling issues, going 
beyond 10 eV or so is formidable and despite advances in 
computing resources seems poised to remain so. The cur¬ 
rent orbital-free extension alleviates this bottle neck, al¬ 
lowing for accurate results in conjunction with the Kohn- 
Sham method, from zero temperature through the high 
temperature Thomas-Fermi limit. 

As a relevant application we constructed a wide rang¬ 
ing EOS for Si02. The EOS constrained by the QMD 
calculations shows very good agreement with the recent 
shock experiment data. An immediate and critical appli¬ 
cation of the new EOS lies in design and analysis of high- 
pressure shock experiments of different materials where 
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a-quartz has been used extensively as a window mate¬ 
rial, and hence determination of these materials relies on 
a highly accurate Si02 EOS. 
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